7  S5: Índices Espectrales

Índices Espectrales con R y Rgee

Author

Centro de Inteligencia Territorial

Published

May 9, 2023

7.1 Introducción

Los índices espectrales son herramientas utilizadas en percepción remota para analizar la reflectancia de la superficie terrestre en diferentes longitudes de onda del espectro electromagnético. Estos índices se basan en la relación entre la reflectancia en dos o más bandas espectrales y se utilizan para medir diferentes propiedades de la vegetación, como la salud, la humedad y la densidad.

Algunos de los índices espectrales más utilizados son:

  • NDVI (Índice de Vegetación de Diferencia Normalizada): Este índice se utiliza para medir la salud y densidad de la vegetación. Se calcula utilizando la reflectancia en la banda roja y la banda infrarroja cercana. Un valor alto de NDVI indica una vegetación densa y saludable.
  • EVI (Índice de Vegetación Mejorado): Este índice también se utiliza para medir la salud y densidad de la vegetación, pero tiene en cuenta factores como la influencia del suelo y la atmósfera. Se calcula utilizando la reflectancia en la banda azul, roja e infrarroja cercana.
  • SAVI (Índice de Vegetación Ajustado para Suelo): Este índice se utiliza para reducir la influencia del suelo en la medición de la vegetación. Se calcula utilizando la reflectancia en la banda roja y la infrarroja cercana y tiene en cuenta la cantidad de luz reflejada por el suelo.
  • NLI (Índice de Línea de Agua Normalizado): Este índice se utiliza para medir la cantidad de agua en una superficie. Se calcula utilizando la reflectancia en la banda verde y la banda infrarroja cercana.
  • NDWI (Índice de Agua Normalizado): Este índice también se utiliza para medir la cantidad de agua en una superficie, pero tiene en cuenta la reflectancia en la banda del infrarrojo de onda corta y la banda verde.

Estos índices espectrales son solo algunos ejemplos de las muchas herramientas disponibles para analizar la vegetación y el agua utilizando percepción remota.

7.2 Índices Urbanos

Algunos ejemplos de índices espectrales normalizados utilizados para analizar el entorno urbano con percepción remota son:

  • NDBI (Índice de Brillo de la Edificación Normalizado): Este índice se utiliza para identificar y mapear áreas urbanas construidas. Se calcula utilizando la reflectancia en la banda del infrarrojo cercano y la banda del rojo.
  • SI (Índice de Suelo): Este índice se utiliza para detectar y mapear la presencia de suelo desnudo en áreas urbanas. Se calcula utilizando la reflectancia en la banda del rojo y la banda del infrarrojo cercano.

Estos índices espectrales normalizados pueden ser útiles para el análisis de la planificación urbana, la gestión del medio ambiente y otros fines relacionados con la ciudad.

Haciendo uso de proyecto awesome spectral indices que ofrece a los usuarios las herramientas para consultar índices espectrales y las bandas requeridas para el cálculo de un índice. Incluso cuentan con Aplicación web

7.3 Calculo de Índices de Urbanos

Cargar Librerías

library(sf)
library(rgee)
ee_Initialize('denis.berroeta@gmail.com', drive = TRUE)
library(tidyverse)
library(ggplot2)
library(viridis)
library(mapview)

Si seleccionamos los índices espectrales obtendremos la siguente lista


Seleccionares índices espectrales urbanos (comparaciones de indices urbanos) y escogeremos el NDBI (Normalized Difference Built-Up Index).

NDBI Normalized Difference Built-Up Index

7.3.1 Selección de Área de Estudio

mi_comuna <-  "LAS CONDES"
roi <-  readRDS("data/rds/zonas_urb_consolidadas.rds") %>% 
  st_as_sf() %>% 
  filter( NOM_COMUNA == mi_comuna) %>% 
  st_transform(4326) 


plot(roi$geometry)

roi_ee <-  roi%>%
  st_union() %>% 
  st_geometry() %>% 
  sf_as_ee()

7.3.2 Selección Imagen Satelital

Harmonized Sentinel-2 MSI: MultiSpectral Instrument, Level-2A

maskS2clouds <- function(image) {
  qa <- image$select("QA60")
  # Bits 10 and 11 are clouds and cirrus, respectively.
  cloudBitMask <- bitwShiftL(1, 10)
  cirrusBitMask <- bitwShiftL(1, 11)

  # Both flags should be set to zero, indicating clear conditions.
  mask <- qa$bitwiseAnd(cloudBitMask)$eq(0)$And(
    qa$bitwiseAnd(cirrusBitMask)$eq(0)
  )

  # Return the masked and scaled data, without the QA bands.
  image$updateMask(mask)$
    divide(10000)$
    select("B.*")$
    copyProperties(image, list("system:time_start"))
}
# Map the function over one year of data and take the median.
# Load Sentinel-2 TOA reflectance data.
S2 <- ee$ImageCollection("COPERNICUS/S2_HARMONIZED")$
  filterDate("2016-01-01", "2016-12-31")$
  filter(ee$Filter$lt("CLOUDY_PIXEL_PERCENTAGE", 20))$
  map(maskS2clouds)$
  median()$
  clip(roi_ee)
viz <- list(bands = c("B4", "B3", "B2"), min = 0, max = 0.3)
Map$centerObject(eeObject = S2,zoom = 12) 
Map$addLayer(S2, viz, "RGB")

7.3.3 Aplicar Índice

la función NDBI corresponde a la siguiente:

NDBI = \frac{S1-N}{S1+N} Por lo anterior es muy importante conocer los nombres de la bandas para crear correctamente la función, a continuación de exponen las bandas del producto satelital sentinel.

F_NDBI<- function (image){
  ndbi<-image$expression('float(S1-N)/float(S1+N)',
                         opt_map=list("N"=image$select("B8"),
                                      "S1"=image$select("B11")))$
    rename("NDBI")
  return(ndbi)
}
S2_NDBI <- F_NDBI(S2)

7.3.4 Visualizar Resultados

pal_ndbi <- c('white' , 'white', 'white', 'red' , 'red', 'red')

viz_ndbi<- list(min=-1,
                max=1,
                palette=pal_ndbi)

Map$centerObject(eeObject = roi_ee,zoom = 12) 
Map$addLayer(S2_NDBI,visParams =viz_ndbi)

7.3.5 Calcular un índice

mz_sii <- readRDS("data/rds/mz_constru_SII.rds") %>% 
  st_as_sf() %>% 
  filter(n_com == mi_comuna) %>% 
  mutate(id = 1:nrow(.))

Extraer los valores con la función ee_extract cuya documentación se encuentra en el siguiente link

#Extract values - getInfo
ee_com_urb <- ee_extract(
 x = S2_NDBI,
 y = mz_sii["id"],
 scale = 10,
 fun = ee$Reducer$mean(),
 sf = TRUE
)
mz_sii_urb <- mz_sii %>% 
  left_join(ee_com_urb %>% st_drop_geometry(), by = "id")

7.3.6 Corrección

mz_sii_urb <- mz_sii_urb %>% 
  mutate(NDBI = ifelse(total == 0, NA, NDBI))

7.3.7 Histograma de NDVI

ggplot(mz_sii_urb, aes(x = NDBI)) +
  geom_histogram(bins = 50, fill = "#1d91c0", color ="gray80")+
  theme_bw()+
  labs(title="Histograma de Índice de Urbanismo", x ="NDBI", y = "Frecuencias")+
  theme(plot.title = element_text(face = "bold",colour= "gray60", size=10)) 

7.3.8 Visualización espacial de NDBI

mapview(mz_sii_urb, zcol= "NDBI")
library(ggplot2)
library(viridis)

# Visualización ggplot y sf
ggplot() +
  geom_sf(data = mz_sii_urb, aes(fill = NDBI), alpha=0.7)+
  scale_fill_viridis_c(option = "E", direction = -1)+
  ggtitle("Urban Index: NDBI") +
  theme_bw() +
  theme(legend.position="none")+
  theme(panel.grid.major = element_line(colour = "gray80"),
        panel.grid.minor = element_line(colour = "gray80"))

7.3.9 Guardar los resultados

saveRDS(ee_com_urb, file = "data/rds/ndbi_2016.rds")

7.4 Monitoreo de Sequía Laguna de Aculeo

Para el siguiente análisis utilizaremos el NDWI: Normalized Difference Water Index)

7.4.1 El Problema

Noticias Relacionadas:

“Laguna de Aculeo, una de las primeras víctimas del cambio climático y megasequía según estudio”. Biobiochile.cl 20-JUN-2019

“La sequía de Aculeo vista desde el espacio”. Diario La Tercera 11-MAY-2018

“En desaparecida laguna de Aculeo proponen que nueva Constitución tenga sello medioambiental: Greenpeace ilustra la grave emergencia climática en Chile” Greenpeace 28-NOV-2019

7.4.2 Lectura de Insumo Cuerpos de Agua Nacional

masas_agua <- st_read("data/shape/masas_agua.shp", quiet = T)
masas_agua
Simple feature collection with 16644 features and 7 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: -5659.382 ymin: 3803587 xmax: 692350 ymax: 8050958
Projected CRS: WGS 84 / UTM zone 19S
First 10 features:
                                    NOMBRE TIPO_MAGUA COMUNA    AREA_KM2 RULEID
1                       Laguna Tebinquiche     Laguna   <NA>  2.08693247      5
2                                     <NA>     Laguna   <NA>  0.84982790      5
3                           Salar de Pujsa      Salar   <NA>  0.22226662      7
4                           Salar de Pujsa      Salar   <NA> 16.64594244      7
5                         Salar de Atacama      Salar   <NA> 49.61221690      7
6                           Salar de Pujsa      Salar   <NA>  0.02636444      7
7            Salar de Loyoques o Quisquiro      Salar   <NA> 65.49493327      7
8            Salar de Loyoques o Quisquiro      Salar   <NA>  0.63644089      7
9                           Salar de Pujsa      Salar   <NA>  0.08644732      7
10 Laguna en Salar de Loyoques o Quisquiro     Laguna   <NA>  0.02395804      5
   SHAPE_Leng  SHAPE_Area                       geometry
1  12353.4677  2086932.47 MULTIPOLYGON (((575205.1 74...
2   4397.5639   849827.90 MULTIPOLYGON (((661560.1 74...
3   3305.8444   222266.62 MULTIPOLYGON (((652931.2 74...
4  26208.2874 16645942.44 MULTIPOLYGON (((652560 7434...
5  44044.4984 49612216.90 MULTIPOLYGON (((569340.4 74...
6    606.6370    26364.45 MULTIPOLYGON (((653547.4 74...
7  67250.4947 65494933.27 MULTIPOLYGON (((672122.6 74...
8   5459.5505   636440.89 MULTIPOLYGON (((671993 7433...
9   1147.0949    86447.32 MULTIPOLYGON (((653825.7 74...
10   686.3645    23958.04 MULTIPOLYGON (((673224.5 74...
#Buscar nombre del cuerpo de aguar por patrón de texto
grep(pattern = "Aculeo$", x = masas_agua$NOMBRE, value = T)
[1] "Laguna de Aculeo"
# Leer Laguna Aculeo
lag_aculeo <- masas_agua%>%
  filter(NOMBRE == "Laguna de Aculeo")%>%
  sfheaders::sf_remove_holes() ## remover poligonos internos%>%>%
plot(lag_aculeo$geometry,col = "skyblue",  main = "Laguna Aculeo")

#Laguna Formato ee
lag_aculeo_ee <- lag_aculeo%>%
  st_buffer(500)%>%
  st_transform(4326)%>%  #crs atlon
  sf_as_ee()

# regiòn  para definir regiòn de interès para los mapas
region <- lag_aculeo_ee$geometry()$bounds()

7.4.3 Lectura de insumo satelital visualización previa

imagen <- ee$ImageCollection('LANDSAT/LC08/C01/T1_TOA')$
  filterBounds(region)$
  filterDate('2018-01-01','2018-12-31')$
  filterMetadata('CLOUD_COVER','less_than', 5)$
  median()$
  clip(region)


ndwi <- imagen$normalizedDifference(c("B3", "B5"))$clip(region)
ndwiMask <- ndwi$updateMask(ndwi$gte(0))#greater than or equal to the given value.



ndwiViz <- list(min = 0, max = 1, palette = c("white", "#00FFFF",'#0080FF', "#0000FF"))
imgViz <- list(min = 0, max = 0.5,  bands = c("B4", "B3", "B2"), gamma = c(0.95, 1.1, 1))

Map$centerObject(region, zoom = 13)
mapa <- Map$addLayer(imagen, imgViz,'RGB')+
  Map$addLayer(ndwiMask, ndwiViz,'NDWI')
mapa

7.4.4 Construcción de un Mosaico con NDWI y Falso Color

imageRGB <- imagen$visualize(bands = c('B5', 'B4', 'B3'), max = 0.5, gamma = c(0.95, 1.1, 1))
ndwiRGB <- ndwiMask$visualize(min = 0, max = 1, palette = c('#00FFFF', '#0000FF'))
mosaic <- ee$ImageCollection(c(imageRGB, ndwiRGB))$mosaic() #función crea mo
Map$centerObject(region, zoom = 13)
Map$addLayer(mosaic , list(), 'mosaic')

7.4.5 Creación de Serie Anual desde el 2015 al 2020

imagen <- ee$ImageCollection('LANDSAT/LC08/C01/T1_TOA')$
  filterBounds(region)$
  filterDate('2015-01-01','2020-12-31')$
  filterMetadata('CLOUD_COVER','less_than', 5)
anual <- ee$List$sequence(2015, 2020)


anual_l8 <- function(y) {
  imagen$filter(ee$Filter$calendarRange(y, y, "year"))$
    median()$
    select("B4", "B3", "B2")$
    clip(region)}

l8_year <- anual$map(ee_utils_pyfunc(anual_l8))
l8_year_2015 <- ee$Image(l8_year$get(0))
l8_year_2016 <- ee$Image(l8_year$get(1))
l8_year_2017 <- ee$Image(l8_year$get(2))
l8_year_2018 <- ee$Image(l8_year$get(3))
l8_year_2019 <- ee$Image(l8_year$get(4))
l8_year_2020 <- ee$Image(l8_year$get(5))

visparams <- list(
  bands = c("B4", "B3", "B2"), min = 0,max = 0.6,gamma = 1.4)

Map$centerObject(region, zoom = 13)
map1 <- Map$addLayer(l8_year_2015, visparams, name = "2015")+
  Map$addLayer(l8_year_2016, visparams, name = "2016")+
  Map$addLayer(l8_year_2017, visparams, name = "2017")+
  Map$addLayer(l8_year_2018, visparams, name = "2018")+
  Map$addLayer(l8_year_2019, visparams, name = "2019")+
  Map$addLayer(l8_year_2019, visparams, name = "2020")
map1
# referencias: # https://github.com/r-spatial/rgee/blob/examples/ImageCollection/creating_monthly_imagery.R